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Abstract 

A videotape presentation of flow in a packed bed 
of spheres is provided. The flow experiment con- 
sisted of three principal elements: ( 1 ) an oil tunnel 
76.2 mm by 76.2 mm in cross section, (2) a packed 
bed of spheres in regular and irregular arrays, and 
(3) a flow characterization methodology, either (a) 
full flow field tracking (FFFT) or (b) computational 
fluid dynamic (CFD) simulation. The refraction 
indices of the oil and the test array of spheres were 
closely matched, and the flow was seeded with alu- 
minum oxide particles. Planar laser light provided 
a two-dimensional projection of the flow field, and 
a traverse simulated a three-dimensional image of 
the entire flow field. Light focusing and reflection 
rendered the spheres black, permitting visualiza- 
tion of the planar circular interfaces in both the axial 
and transverse directions. Flows were observed near 
the wall-sphere interface and within the set of 
spheres. 

The CFD model required that a representative 
section of a packed bed be formed and gridded, en- 
closing and cutting six spheres so that symmetry 
conditions could be imposed at all cross-boundaries. 
Simulations had to be made with the flow direction 
at right angles to that used in the experiments, how- 
ever, to take advantage of flow symmetry. Careful 
attention to detail was required for proper gridding. 
The flow field was three-dimensional and complex 
to describe, yet the most prominent finding was flow 
threads, as computed in the representative ‘cube’ 
of spheres with face symmetry and conclusively 
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demonstrated experimentally herein. Random pack- 
ing and bed voids tended to disrupt the laminar flow, 
creating vortices. 

Nomenclature 

A area 

D sphere diameter 

Go mass flux (no spheres in tunnel) 

L length, model of spheres 

Ap pressure drop 

Re Reynolds number 

V volume 

x,y,z coordinates 

E volume porosity 

p density 

ji viscosity 

Introduction 

Although similar visualization materials and tech- 
niques have been discussed in previous work by 
Hendricks [1,2], it was realized that printed frames, 
grabbed from videotapes of the flows in beds of 
packed spheres failed to reproduce what was visu- 
alized either in computational fluid dynamic (CFD) 
virtual images or experimentally. And herein, the 
essences of the numerical and experimental work 
are synergistically combined to provide background 
for the visualization that can be realized only by 
viewing the videorecording (of the experiment or 
the CFD virtual images) or the experiment itself. 

Packed spherical particle beds are characterized, 
for example, by Ergun [3] and Bird et al. [4] as em- 
pirical combinations of laminar and turbulent flows. 
Process design engineers use similar methods to 
define packed beds of spheres for a wide variety of 
porous media applications (filters, catalytic packs, 
heat exchangers, etc.). Yet the details of heat, mass, 
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and momentum transfers in these packed structures 
are not understood. Studies made of a dumped bed 
of spheres by Jolls and Hanratty [5] and confirmed 
by Karabelas [6] show that the flow is three- 
dimensional and laminar at Reynolds numbers (Re) 
less than 40, with unsteady flows occurring at 
1 10 < Re < 150. Flows in cubic arrays of spheres 
are steady to Re = 82 and characterized by nine 
regions of reverse flow with similar patterns 
observed for Re = 200, even though the flow 
becomes unsteady [7]. In numerical modeling of 
porous media Saleh et al. [8], Rosenstein [9], 
and Yarlagadda and Yoganathan [10] found that 
little knowledge of the interstitial flow fields was 
available. Saleh used particle image planar displace- 
ment velocimetry and integration of the continuity 
equation to describe the third component (similar 
to the procedure used by Rosenstein). 

The full flow field tracking (FFFT) method 
[11,12] can be applied to visualize and quantize 
the flow patterns and fluid velocities within a 
packed bed of spheres. The numerical simulations 
were made with the code SCISEAL, a pressure- 
based, finite-volume, three-dimensional CFD code 
with implicit multiple-block capability [13]. More 


details of packed-bed sphere flows can be found in 
Athavale et al. [14]. The gridding was executed by 
using CFB-GEOM [15]. Materials from Athavale 
et al. [13,14] and Hendricks et al. [1,2] are used to 
orient the reader in viewing the videotape. 

Test Facility 

The test facility consisted of an oil tunnel, flow sys- 
tem components, videorecording equipment, a 
4-W continuous-wave (4CW) argon-ion laser, lens 
systems, data recorders, and a test configuration 
(fig. 1). The matrix of spheres (fig. 2) and entry 
ramp spacer was fabricated to fill the tunnel cross 
section (76.2 mm by 76.2 mm). Both the laser light 
sheet and the viewing port walls were Lucite. The 
test sphere matrix was fabricated by using optical 
quality glue to bind the array, matching the refrac- 
tion index of Lucite, 1.4905. This technique sacri- 
ficed the point of contact, yet permitted construc- 
tion of a regular array of staggered spheres and a 
clear view through the spheres when they were 
immersed in the oil, which also had the same re- 
fraction index as Lucite. The spheres in the bottom 
layer were ground flat, and their height was ad- 
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Figure 1 . — Schematic of test facility. 
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x = 0 12 3 4 5 6 7 8 


Figure 2. — Matrices of Lucite spheres of refractive index 1 .4905; sphere diameter, 1 2.7 mm. x = 0(front 
wall), x = 1 (6.4 mm), x = 2 (12.7 mm), x = 3 (19.1 mm), x = 4 (midplane), (a) Side view (camera view), 
(b) Front view (flow view). 


justed to exactly fit the tunnel. No other spheres 
were altered. The width was adjusted by using a 
gradual entry ramp that blocked a portion of the 
tunnel and at the same time eliminated any recir- 
culating flow at the entry of the array. 

The packed sphere bed (78.7 mm long by 
50.8 mm wide by 76.2 mm high) contained 
172 spheres. The 7-by-7 longitudinal array 
(fig. 2(a)) had a transverse array of four 12.7-mm- 
diameter spheres (fig. 2(b)). The volume of the test 
array envelope was V ( = 305 cm 3 , the solid volume 
of spheres was V s = 1 84 cm 3 . The void V Q =V t - V s 
= 121 cm 3 , with an effective open cross section for 
flow of A e = 17.9 cm 2 . The geometric arrangements 
for the numerical analysis and sphere identifica- 
tion are given as figures 6 and 7 and discussed later. 

The coherent beam of the 4CW argon-ion laser 
was directed by micrometrically adjustable mirrors 
through two cylindrical lenses positioned at 90° to 
each other, through the Lucite tunnel window, and 
into the test section. The light sheet was approxi- 
mately 0.1 mm thick, and the flow was seeded with 


magnesium oxide flow tracers. Micrometric adjust- 
ments controlled the position of the light slices 
scanning across the test section and provided, 
through composition, a three-dimensional visual- 
ization of the flow field. The light was naturally 
focused and reflected at the boundary interfaces 
between the working fluid and the spheres, render- 
ing the spheres black and permitting visualization 
of the planar circular interfaces in both the axial 
and transverse directions. The images were 
recorded on a videocassette recorder. 

Other instrumentation included static pressure 
taps and thermocouples upstream and downstream 
of the test matrix and two spheres coated with liquid 
crystal with a narrow band response about 21 °C. 

Experimental Results 

Figure 3 represents the flow field 6.35 mm from 
the front window (x = 1) in the array of 12.7-mm- 
diameter spheres. The flow rate was 3.8 1/min 
(1 gal/min). The spheres were numbered to facili- 
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Figure 3. — Visualization of flow in 1 2.7-mm-diameter sphere matrix (fig. 2) 6.4 mm from front wall (x = 1 .0) with 
flow at 3.8 liters/min (1 gal/mi n). 


tate flow descriptions. The laser light slice was 
strobed at 180 Hz, giving a six-segment particle 
streak per recorded videotape frame. The inlet flow 
streaks around spheres 1 to 4 resemble flows around 
an array of cylindrical pins. Stagnation zones 
occurred in front of and behind spheres 1 to 4, and 
wake zones appeared behind spheres 22 to 25. Flow 
near the tunnel upper wall, spheres 4, 11, 18 and 
25, showed a high level of distortion and lower 
velocities, illustrating the influence of a sphere- wall 
boundary. Careful attention to the particle trace 
lengths and brightness defined a flow thread in the 
axial direction, which for regular arrays resembled 


a complex sine wave. The flow spiraled upward 
from the gap between spheres 13 and 16, then 
slowed arid turned downward into the pocket 
formed by spheres 16, 19, 20, and 23, and flowed 
out betwee n spheres 20 and 23. The streak between 
spheres 1 and 2 headed to the stagnation zone of 
sphere 5. The wake of sphere 5 was vivid and a 
streak extended toward the stagnation zone of 
sphere 12. The wake behind sphere 14 and the stag- 
nation zone in front of sphere 21 were almost con- 
nected. Flew symmetry existed between spheres 10, 
13, 14, and 17 and the closing of small wakes behind 
spheres 2 and 3. The flow thread was complex. 
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Figure 4.— Transverse visualization of flow in 1 2.7-mm-diameter sphere matrix (fig. 2) at column of spheres 5, 6, 
and 7 (fig. 3) with flow at 1 1 .4 liters/min (3 gal/min). (a) x = 0 (front wall), (b) x = 1 (0.5 sphere diameter from 
wall), (c) x = 2 (1 sphere diameter from wall), (d) x = 3 (1 .5 sphere diameters from wall). 

Midway across the matrix block, 25.4 mm from weaving flow thread are partially visible. In fig- 
the front wall (x = 4), the effects of the upper wall me 4(c), at x = 2 (12.7 mm, or 1 sphere diameter), 

on the flow between spheres 1 1, 14, 18, 21, and 25 the light sheet is skewed and portions of the adja- 

was still pronounced, with diminished velocities cent spheres and their multiple contact points 
between spheres 18, 21, and 25 and the wall. become visible. Figure 4(d), atx = 3 (19.1 mm, or 

At 1 1 .4 1/min (3 gal/min) the boundary layers, 1 .5 sphere diameters), shows the flow thread around 

jetting, and wakes became more pronounced. Fig- a sphere and adjacent contact points, 
ure 4 is a transverse scan of flow that includes the Figure 5 presents a closeup view of the flow in 

column of spheres 5, 6, and 7 (see figs. 2 and 3). the upper part of figure 4(d) at a magnification of 

Figure 4(a) shows low- velocity flow streaks and 7.8. The strong effects of the flow are apparent from 

threads about the sphere in the immediate vicinity a stagnation region in the cusp at the confluence 
of the front lateral wall (x = 0). Figure 4(b) shows between the spheres and the extremely slow flow 
higher velocity flow streaks at x = 1 .0 (6.4 mm, or in the channel between the spheres. 

0.5 sphere diameter), and particles tracking the 
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Figure 5. — Sphere matrix flow image at 7.8X for 12.7-mm-diameter sphere matrix 
and 1 1 .4-liters/min (3-gal/min) flow rate. 



Ce 


Figure 6. — Stacking of sphere layers to generate densely packed arrangements. 
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Numerical Modeling 

Unlike visualization of flows in pinned arrays, 
which are weakly three-dimensional dependent, the 
flow, heat, and mass transfers in spherical beds are 
strongly three-dimensional dependent and require 
the combined efforts of numerical and experimen- 
tal techniques to embody and characterize models 
and results. Also in contrast to pinned arrays, the 
geometric descriptions become quite complex. 
Spherical clusters with hexagonal or cubic densely 
packed beds have 12 neighbors and a void fraction 
of about 0.26 (figs. 6 and 7). The first-ever simula- 
tions of the flow in the interstitial spaces of a hex- 
agonal densely packed arrangement were done by 
Athavale et al. [14] using structured, multiple- 
domain grids. T\vo types of interstitial space were 
identified as the basic building blocks for generat- 
ing the flow by replication. The two types were (1) 
tetra-space, bounded by four spheres with centers 
at the apex of a tetrahedron; and (2) hexa-space, 
bounded by six spheres with centers at the apexes 
of a double pyramidal structure. Each of these 
spaces was further subdivided into simpler domains 
so as to be able to build structured grids in the flow 
domain. By joining the tetra- and hexa-spaces at 
appropriate planes, a given volume could be popu- 
lated with the closely packed spheres. Flow and 
heat transfer simulations on these spheres were re- 
ported by Athavale et al. [14], The complexity of 
the flow domain, however, limited the overall do- 
main size that could be considered with structured 
grids. The number of domains and grid cells grew 
rapidly with each tetra- and hexa-space added to the 
flow problem. Additionally, each space contained 
a singularity point, and the grid topology often forced 
grid clustering in areas where it was not needed, 
leading to nonoptimal use of the available cells. 

With an unstructured grid format, where the 
basic volume element is a tetrahedron, the com- 
plex shapes of the interstitial spaces could be 
gridded more easily, with better control of surface 
and volume discretization in terms of cell sizes and 
clustering. In addition, with the unstructured grid 


'For identifying the basic spaces ping-pong balls are useful 
mnemonics. 


format the representative ‘slice’ of the flow domain 
could be easily created by using much simpler cut- 
ting planes (e.g., parallel to xy, yz, and zx planes in 
the Cartesian coordinate system) rather than the 
naturally occurring interfaces that had to be used 
with the structured grids. Finally, by using proper 
cutting planes and locations, symmetry and 
periodicity boundaries could be used to make the 
flow in the representative slice close to the real situ- 
ation in a packed bed. 

Figure 8(a) 1 shows a solid model of such a rep- 
resentative slice, with surface triangulation on all 
outer boundary surfaces (fig. 8(b)). The curved sur- 
faces are all outer surfaces of the bounding spheres 
and are no-slip walls. The planar surfaces in be- 
tween are the areas where flow can take place. With 
this grid orientation, for a direct comparison with 
the experiments, the bulk flow should be along the 
z direction. However, because the x z plane surfaces 
are not in appropriate geometric locations for sym- 
metry conditions, the flow direction was oriented 
along the y direction to allow the use of symmetry 
conditions on the cross-plane surfaces (yz and yx 
planes), which are identified in figure 8(a). With 
this change the viewing plane still remained the 
same (yz plane), but the flow direction in the simu- 
lations was at right angles to the experimental 
direction. Differences in the numerical and experi- 
mental cube of spheres are shown sketched in fig- 
ures 9(a) and 9(b) respectively with tabulated 
porosities listed in table I. Although this prevented 
a direct comparison with the experiments, the simu- 
lations still showed a great deal of similarity with 
the experiments and offered insight into the com- 
plexity of the flow field. To allow bulk flow in the 
z direction, periodic conditions would need to be 
imposed in all three directions, and this work is 
under way. 

The main flow direction, then, is along the y 
coordinate axis as shown. A uniform flow velocity 
was prescribed on the ‘open’ areas of the inlet xz 


Table I. — Porosities of numerical and 
experimental cube of spheres models. 


Volume porosity. 
Face porosity 

Numerical 

0.260 

Experimental 

0.302 

xy 

0.093 

0.093 

xz 

.198 

.093 

^ 

.444 

.476 
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pressure [2,3,4] 
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Sphere 4 

Specified pressure 
[1,3,4] 
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Inlet 
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s 

Uniform inlet 
flow conditions 


Sphere 2 surface grid 
shown on solid wall to 
show details 


Symmetry 

[1,2,3] 


Nomenclature 
[1 ,4,6] * plane cutting 

spheres 1 ,4, and 8 


zi 


y 
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parallel 
to x-y 
plane 


Nomenclature 
[1 ,4,6] = plane cutting 

spheres 1 ,4, and 6 


Sphere 6 


Sphere 4 
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shown for 
clarity 


Inlet 

flow 


(b) 


' S-Jj 

[1,4,6] 

[1,2,4] 

[2,4,5] 

f Uniform inlet 
flow conditions 

Sphere 1 


[1,3,6] 

Specified pressures 
' - [3,5,6] behind sphere 5 

Sphere 5 


Sphere 3 



Sphere 2 


Specified 
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[2,3,5] 

behind 

spheres 

3 and 2 


Specified pressures 
[1 ,2,3] behind sphere 
2 parallel to x-y plane 


Figure 7. — Solid model, showing flow geometry and boundary conditions, (a) Single tetra-space. (b) Single 
hexa-space. 
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Figure 9.— ‘Cube’ of spheres sketches, (a) Numerical, 
(b) Experimental. 


plane, and a constant static pressure was specified 
at the opposite xz plane, the exit boundary. The 
uniform velocity and pressure values specified on 
the xz boundaries introduce ‘entry’ effects into the 
flow field for the flow Reynolds numbers consid- 


2 It has been shown that for an absolutely regular, low- 

porosity (tightly packed) array of cylinders turbulence can 

not be engendered [16]. 


ered. However, these effects were expected to last 
over a small distance inward from these two bound- 
aries, and the flow in the central portion of the do- 
main can be compared with that in the experiments. 

Incompressible, constant-property, laminar 
flow 2 of oil through the matrix was simulated at the 
two volume flow rates used in the experiments. The 
inlet velocities corresponded to flow rates of 3.8 
and 1 1 .4 1/min ( 1 and 3 gal/min) in the test appara- 
tus described above. The sphere diameters were set 
to 12.7 and 19.05 mm. The oil had a densityof 900 
kg/m 3 and a kinematic viscosity of 1 .0x1 O' 4 m 2 /s. 
The free tunnel inlet velocities (no spheres in the 
tunnel) were 0.016 and 0.048 m/s, respectively, and 
adjusted by the flow cross-sectional area (with 
spheres in the tunnel) to yield the flow rates used 
in the analysis; A total /A open | y=0 = 5.04. 

These and other flow rates are discussed in 
Numerical Results. 

The flow solver used for simulations was 
SCISEAL(U), which has unstructured/hybrid grid 
capability. The code uses a pressure-based, finite- 
volume integration method to solve Navier-Stokes 
equations. It can be used to treat compressible and 
incompressible flows and has a variety of bound- 
ary condition types and conjugate heat transfer and 
turbulence models [13]. 

Numerical Results 

Figures 10(a) and (b) are representative velocity 
vector plot s within the array of spheres at two con- 
stant x planes for a flow rate of 3.8 1/min (1 gal/ 
min). The orientation and magnitude of the veloc- 
ity formed the flow thread as the yz cutting plane 
was moved one sphere diameter across the flow 
field (experimental flow is along the vertical direc- 
tion in the se figures). Although flow orientation, 
grid refinement, and postprocessing of the results 
require further resolution, initial values of the 
numerical pressure drop across the ‘dual cube’ were 
comparable to those calculated from the Ergun 
relation [2,4] (table II and fig. 11). Approximate 
average values of the static pressures over the inlet 
and outlet planes were used to calculate the pres- 
sure drop over the ‘dual cube’ and then normalized 
with the length of the flow domain in the flow 
direction to get the linear pressure drop rates. 
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(a) 



(b) 

Figure 10. — Velocity vectors in single-axis hexa-space along yz (x = constant) cutting 
plane in sphere cluster with flow rate at 3.8 liters/min (1 gal/min). (a) x = 4. (b) x = 5. 


Table II. — Comparison of numerical and calculated Ergun [3,4] pressure 
drops across arrays of spheres (Ergun parameters: porosity = 0.26, 
density = 900 kg/m 3 , viscosity = 0.0904 Pa-s). 



Inlet velocity, m/s (flow rate, g 

al/min) 


0.016 

(1) 

0.048 

(3) 

0.08 

(5) 

0.144 

(9) 

0.288 

_Q8) 

0.576 

(36) 

1.152 

(72) 

| Pressure drop, MPa/m 1 

Sphere diamter, 
12.7 mm (1/2 in.) 
Ergun relation 

0.042 

.043 

0.138 

.138 

0.249 

.243 

0.507 

.486 

1.25 

1.19 

3.34 

3.24 

9.94 

9.95 

Difference, percent 

3 

0 

-2.2 

-4.3 

-5.1 

-2.9 

0.1 

Sphere diameter, 
19.05 mm (3/4 in.) 
Ergun relation 

0.0192 

.01% 

0.065 

.064 

0.12 

.116 

0.252 

.24 

0.653 

.625 

1.848 

1.827 

5.738 

5.964 

Difference, percent 

1.8 

-1.5 

-3.8 

-5.1 

-4.5 

-1.2 

3.8 
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Data Numercai Prior 

Braun et al. [14] results results 



Figure 11 . — Numerical and experimental results superimposed on a sketch illustrating the general behavior of 
the Erg un-relation (Bird et al. [4]). 


Departures of the experimental data from the Ergun- 
locus are most likely due to difficulties in measur- 
ing the local pressure drops at the appropriate 
positions within the array and require further in- 
vestigations. These figures, the video recording, and 
the experiments reveal that clear flow threads were 
established by the interstitial flow domain shapes. 
In general, the flow patterns were similar for the 
two flow rates and correlated well with the experi- 
mental pattern. The flow velocities (Reynolds num- 
bers) were small enough that there were no major 
recirculation zones at tunnel inlet velocities below 
0.6 m/s (flow rates below 36 gpm); incipient 
computational instabilities and vortical structures, 
figure 12, occurred above 1.15 m/s (flow rates above 
70 gpm). These transition effects are being investi- 
gated. Some evidence of secondary vortical struc- 
tures in the flow field was seen in the earlier flow 


simulations by Athavale et al. [14] at fairly high 
Reynolds numbers (over 400). 

Thermal Simulation Using Thermal Liquid 
Crystal (TLC) 

For a set of experiments two liquid-crystal-coated 
spheres were introduced into an ordered but not 
regularly packed bed. Dissipated pump power 
slowly healed the circulated fluid, and distinct color 
banding illustrated regions of high thermal 
gradients rear points of contact. The flows in these 
regions varied from zero to the free-stream flow 
thread, but near the surface the velocities were low. 
Upstream, the heat transfer was by convection; near 
the points of contact, conduction dominated; and 
downstream, wakes and rapid changes in curvature 
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Recirculation zones 

Figure 12a. — Velocity vector plots 12.7-mm-diameter ‘cube’, flow rate at: 
(i) 136 liters/min (36 gal/min). (ii) 273 liters/min (72 gal/min). 


provided mixed heat transfer. Distinct thermal 
variations outlined these regions. 

Videorecording 

The complexity of the flow field, whether virtual 
or experimental, became vivid through the video- 
recording. Visualization of the flow fields revealed 
flow threads, wakes, stagnation zones, and the 
strong influence of the wall-sphere interfaces. The 
tape is available in 1/2-in. VHS format. The 
essences of both the theoretically predicted flow 
threads and those observed experimentally can be 
visualized during a scan of the flow field from front 
to rear lateral walls. These flow threads weave 
through the packed array of spheres in the lateral 
direction, the direction of the bulk flow. Details of 
the flow boundary layer close to the wall and pro- 


gressing through the matrix to the opposite wall 
are revealed. The videorecording can also be used 
to determine quantitative experimental information, 
such as the flow velocities, by using the FFFT tech- 
nique [11,12] and comparing these to the virtual 
flows determined by the CFD solutions. 

Summary 

A videotape presentation is provided that enables 
the viewer to compare experimental and virtual 
computational fluid dynamic flows through packed 
beds of spheres or porous media. For the experi- 
ment the refraction indices of the oil working fluid 
and the spheres were closely matched so that only 
the boundaries of the physical spheres remained 
faintly visible to the eye and the cameras. The sys- 
tem fluid dynamics were made visible by seeding 
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Figure 12b. — Velocity vector plots 19. 05- mm -diameter 'cube', flow rate at: 
(i) 136 liters/min (36 gal/min). pi) 273 liters/min (72 gal/min). 


the flow with aluminum oxide particles and illu- 
minating the packed bed with a planar laser light 
sheet aligned in the direction of the bulk flow. The 
two-dimensional projection of the flow field was 
recorded at right angles to the bulk flow. As the 
planar light sheet traverses from the front wall to 
the rear wall of the tunnel and returns, completing 
one cycle, a three-dimensional image of the entire 
flow field can be visualized. The boundary inter- 
faces between the working fluid and the spheres 
rendered the spheres black, permitting visualiza- 
tion of the planar circular interfaces in both the axial 
and transverse directions and providing direct 
visualization of the flow in complex interstitial 
spaces between the spheres within the bed. For the 
CFD virtual flows, simulations were performed 
with the flow direction at right angles to that used 
in the experiments to allow symmetry plane speci- 
fications. Even with this alteration for a ‘cube’ with 
face symmetries the flow patterns and threads were 


similar to those found experimentally. Simulations 
with periodic conditions and higher tunnel inlet 
velocities are under way that will allow exact 
comparison. 

The pre ssure drops per unit length in the packed 
bed were also calculated and compared with 
results from Ergun’s theory for porous media, and 
good correlation between the two sets was obtained. 

Flows were observed near the plane surfaces 
and within the set of spheres. The computational 
bed was regular, but the experimental bed was not 
as regularly aligned or uniformly ordered. Although 
the flow is truly three-dimensional and complex to 
describe, the most prominent finding was flow 
threads as computed in the representative ‘cube’ of 
spheres with face symmetry and conclusively dem- 
onstrated experimentally herein. The bulk of the 
flow field always tended to find natural flow paths, 
called threads, through the packed bed that are par- 
allel and distinct. For a regular array the number of 
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